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ABSTRACT:  This  paper  is  concerned  with  a  method  of  finite  differences  for 
determining  two-dimensional  and  axisymmetric  supersonic  nozzle  contours.  The 
approach  taken  is  to  specify  a  Mach  number  or  velocity  array  along  the  entire 
centerline  of  the  nozzle  and  then  to  integrate  the  equations  numerically  to 
obtain  the  desired  nozzle  shape.  In  spite  of  the  fact  that  the  original  problem 
is  not  "well  posed"  in  the  subsonic  region,  reasonable  results  were  found 
provided  the  Mach  number  gradient  was  not  too  steep  in  a  neighborhood  of  the  sonic 
line. 
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CHAPTER  I 

INTRODUCTION 

The  study  of  two  dinensional  and  axisymmetric  converging-diverging  nozzles 
has  received  a  great  deal  of  attention  by  both  engineers  and  mathematicians.  Its 
uses  in  wind  tunnels  and  rockets  make  it  of  practical  importance  whereas  its  mixed 
mathematical  nature  give  it  theoretical  interest. 

A  converging-diverging  nozzle  (see  Figure  1)  consists  of  a  converging  section 
where  the  flow  is  subsonic  and  the  mathematical  equations  are  elliptic,  a  diverging 
section  where  the  flow  is  supersonic  and  the  equations  are  hyperbolic,  and  an 
intermediate  transonic  section,  called  the  nozzle  throat,  where  the  flow  passes 
from  subsonic  to  supersonic  and  both  elliptic  and  hyperbolic  behavior  are  present. 
In  wind  tunnel  applications  there  is  a  fourth  section,  the  test  section,  where  the 
flow  is  uniform. 

There  are  basically  two  approaches  taken  when  investigating  nozzles  -  direct 
and  indirect.  In  the  direct  method  a  nozzle  contour  1l  prescribed  and  one  seeks 
to  determine  the  flow  field  inside.  In  the  indirect  method  the  flow  along  the 
axis  of  symmetry,  called  the  centerline,  is  prescribed  and  the  nozzle  giving  rise 
to  this  flow  is  sought.  In  the  design  of  nozzles  the  natural  method  to  use  is 
the  indirect  method  and  thus  will  be  what  we  shall  consider  here.  For  information 
on  direct  methods  the  reader  is  referred  to  [  3  ]  and  [  9  ]. 

The  two  primary  indirect  methods  in  use  are  the  series  methods  and  the  method 
of  characteristics.  In  the  series  methods,  originated  by  Friedrichs  [4  ], 

(see  also  [5]),  one  assumes  that  the  variables  of  Interest  can  be  expanded  in 
series  involving  the  stream  variables.  These  are  then  substituted  into  the  equa¬ 
tions  to  determine  the  values  of  the  coefficients  in  the  expansions.  The  series 
are  then  truncated  and  the  remaining  finite  number  of  terms  is  used  to  compute 
approximations  to  the  solutions.  The  method  of  characteristics  ([10],  [ 11] ) 
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utilizes  a  change  of  variables  to  characteristic  coordinates.  This  transformation 
makes  the  two  independent  variables  in  the  original  equation  additional  dependent 
variables.  One  is  thus  led  to  a  system  of  four  partial  differential  equations  in 
four  unknowns.  The  simplification  occurs  in  that  each  of  the  four  new  equations 
contains  differentiations  with  respect  to  only  one  of  the  new  independent  variables. 
This  enables  one  to  determine  the  characteristics  and  the  flow  in  the  nozzle. 

This  method,  of  course,  is  restricted  to  the  supersonic  region  where  the  equations 
are  hyperbolic  and  the  characteristics  are  real. 

The  great  difficulty  in  the  series  methods  is  the  lack  of  mathematical 
justification.  Little,  if  any,  has  been  accomplished  in  analyzing  the  size  of  the 
terms  which  are  dropped.  The  method  of  characteristics  on  the  other  hand  can  be 
used  only  in  the  supersonic  region.  To  get  around  this  difficulty  different 
methods  e  used  in  each  of  the  different  sections  and  the  solutions  obtained 
are  then  "patched"  together  with  a  "French  curve."  Thus,  for  example,  one  might 
use  the  method  of  characteristics  in  the  supersonic  region  and  patch  the  solution 
so  obtained  to  the  solution  obtained  in  the  subsonic  and  transonic  regions  by 
the  series  method  or  by  a  quasi-one-dimensional  flow  approximation.  However,  in 
many  applications  such  as  very  high  Mach  number  wind  tunnels  where  heat  transfer 
in  the  throat  of  the  tunnel  is  sizable  and  in  short  nozzles  where  the  two 
dimensional  effects  cannot  be  ignored,  the  method  of  patching  leaves  much  to  be 
desired.  It  is  thus  desirable  to  investigate  furthet  methods  for  the  determination 
of  nozzle  contours  for  prescribed  centerline  conditions. 

In  this  paper  we  consider  an  indirect  method  which  is  usable  from  the 
subsonic  region  to  the  supersonic  test  section  of  constant  state.  The  procedure 
is  to  prescribe  a  Mach  number  or  velocity  distribution  along  the  axis  of  symmetry 
and  then  to  integrate  the  governing  equations  of  motion  using  finite  differences. 

We  do  not  claim  to  present  here  a  mathematically  rigorous  technique  but  rather 
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offer  another  method  of  calculation  which  has  advantages  and  disadvantages  when 
compared  with  other  methods.  In  fact  it  will  be  shown  in  the  next  chapter  when 
we  derive  the  equations  of  motion  that  this  procedure  is  unstable  in  the 
subsonic  region  (as  one  would  expect) .  We  shall  say  more  about  this  in  Chapter  IV 
where  we  discuss  our  results  and  conclusions.  The  primary  advantage  of  this  type 
of  method  is  that  it  provides  a  uniform  scheme  from  the  subsonic  region  to  the 
test  section,  thus  eliminating  the  need  for  multiple  patching.  The  disadvantages 
will  be  discussed  in  Chapter  IV. 
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CHAPTER  II 

DEVELOPMENT  AND  STABILITY  OF  THE  DIFFERENTIAL  EQUATIONS 
We  suppose  in  the  following  that  we  are  dealing  with  a  steady,  irrotational, 
perfect  gas  with  constant  specific  heats.  We  consider  both  two-dimensional  and 
axisymmetric  flows.  For  such  flows  the  equation  of  continuity  takes  the  form 

f3(pu)  ,  3(pv)  .  kpv*|  _  n 

3x  3y  y  J 


(2.1)  7  ■  (£\-  — 


where  k  ■  0  for  two-dimensional  flow  and  k  ■  1  for  axisymmetric  flow.  The  equation 
of  state  is 


and  Bernoulli's  equation  can  be  written 


An  explanation  of  the  notation  can  be  found  on  page  18 ,  and  a  derivation  of  the 
above  equations  from  the  basic  equations  of  fluid  dynamics  can  be  found  in  any 
standard  book  on  gasdynamics  (see  e.g.,  [12]) • 

Equation  (2.1)  implies  the  existence  of  a  stream  function  ip  such  that 


(2.4) 


3y 


P_ 


uy 


k  m  p 


3x 


vy 
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and  che  assumption  of  irrotationality  implies  the  existence  of  a  velocity  potential 
$  such  that 

-»■ 

q  ■ 


or 

(2.5)  u-|i  ,  v-|i 

3x  ’  3y 

Substituting  (2.5)  into  (2.4)  we  obtain  the  equations 

(2.6)  ykii  _  _Lk  1± 

3y  p  y  3x  *  3x  p  7  3y 
0  0 

Finally,  from  (2.2)  and  (2.3)  it  follows  that 

1 

(2.7)  ^-(1+^mV”1 

P  2 

We  next  introduce  a  set  of  dimensionless  variables  according  to  the  following 
definitions : 


(2.8)  x* 


x 

L 


r 


-  JL 

a  L 

o 


•  ■ 


a  L 
o 


k+1 


P' 


- 

po 


v,'.-.are  L  is  some  standard  length.  The  previous  equations  then  become  (with  the 
primes  dropped  for  convenience) 
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(2.9) 


11 . 

3y 


k  34 

py  st 


p  ■  (1  + 


3jj) 

3x 


izi 


k  34 

-py  37 


M2) 


1 

Y-l 


Interchanging  the  role  of  dependent  and  Independent  variables  (which  Is 
permissible  since  the  Jacobian  of  the  transformation, 

(2.10)  J(^)  -  pyk(u2  +  v2), 

is  nonzero  off  the  axis)  and  Introducing  new  variables  £  and  n  defined  by 

(2.11)  $  -  /  q(x)dr  ,  ip  -  “  p*q%k+1 

(2.9)  becomes 


(2.12) 


*  *  k 

lx  .  £_a_  (R)  1* 

3n  pt  y  3C 


*  *  k 

3*  .  P  9  /Ih  IX 
3n  pT  y  H 


p 


(l 


q2u> 

2  J 

xc  +  y( 


1 

Y-l 


*  * 

where  (f  denotes  the  velocity  magnitude  along  the  axis  of  symmetry  and  p  and  q 
are  the  sonic  density  and  velocity  magnitude  which  are  fixed  for  given  stagnation 
conditions. 

It  is  the  system  (2.12)  which  we  integrate  to  find  the  nozzle  contour.  We 
view  the  problem  as  an  initial  value  problem  with  information  imposed  for 
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4 


n  ■  0,  the  nozzle  centerline,  and  wish  to  determine  the  solution  for  positive 
n.  We  require  along  n  ■  0  that  y  -  0  and  x  ■  5.  We  furthermore  prescribe  a 
Mach  number  distribution  or  equivalently  a  velocity  distribution  along  this  line. 
This  information  determines  x,  y,  and  p  and  the  derivatives  of  x  and  y  on  the 
centerline  using  (2.12).  This  Is  straightforward  except  perhaps  the  determination 
of  -|^  in  the  axlsymmetrlc  case,  k  -  1,  where  we  must  remove  the  indetermlrancy  of 
~  along  the  centerline.  This  is  done  by  replacing  ^  by  to  obtain  the 
relationship 


h. 

3n 


*  * 


As  the  lines  n  "  constant  correspond  to  streamlines,  any  positive  n  serves 
as  a  nozzle  boundary.  The  particular  n  chosen  to  terminate  the  integration  li 
dete ’.mined  by  various  design  or  convenience  considerations.  The  numerical 
procedure  used  to  perform  this  Integration  Is  discussed  in  the  next  section. 

We  conclude  this  section  with  an  analysis  of  stability  as  suggested  by 
Von  Neumann  and  Rlchtmyer  [13].  We  suppose  that  a  perturbation  fix,  fiy,  fip 
Is  Introduced  and  attempt  to  see  the  effect  of  this  perturbation  at  points  away 
from  the  centerline.  We  thus  replace  x  by  x  +  fix,  y  by  y  +  fiy,  and  p  by  p  +  fip 
In  (2.12)  and  examine  the  effects  of  such  a  perturbation. 

The  equations  of  first  variation  are 

*  * 

pyk-ji-  (fiy)  +  kpy^fiy  +  ykyn<5p  -  n  (fix)  -  0 

A  A 

(2.13)  Pyk  (5x>  +  xnykfiP  +  kxnpfiy  +  nk  ~  (5y)  -  0 

■2xe(pY_1-l)|^(fix)  +  ly^pY"1-!)  ^  (fiy)  +  (Y-l)pY*2(x2  +  y2)fip  -  0 
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Equations  (2.13)  are  systems  of  differential  equations  for  unknowns  6x,  6y, 
and  6p.  The  coefficients  in  the  differential  equations  depend  upon  x,  y,  and 
P  and  are  considered  to  be  smoothly  varying  with  £  and  n.  As  numerical  errors 
behave  like  rapidly  varying  perturbations  we  view  the  coefficients  in  (2.13)  as 
being  constant  In  a  small  region  and  seek  solutions  of  the  form 

rt  i,\  .  i8£  +  on  .  x  i@ £  +  c»n  ►  .  iB£  +  on 

(2.14)  dx  -  6Xq  e  ,  6y  ■  6yg  e  ,  6p  ■  6pg  e 

where  fix^,  dy^,  Pq>  a,  and  6  are  constants  and  p  is  real  and  large.  Putting 
(2.1<i)  into  (2.13)  results  In 

*  * 

(pyko  +  kpy^)  6yQ  +  yN^Pg  “  ^  ±6p«»5  *  0 

*  * 

(2.15)  (kpx^  +  4^-  n^iB)  6yQ  +  x^SpQ  +  p  yk  a6xQ  -  0 

2y^(pY"1-l)  165y0  +  (y-D  PY~2(x2  +  y2)  <SPq  +  2x?(p1,-1-l)iB6x0  -  0 

We  are  thus  led  to  two  systems  of  linear  homogeneous  equations  In  6xq,  6y^, 
and  5pQ.  Assuming  the  perturbations  are  nontrivial  we  must  have  that  the  coeffi¬ 
cient  determinants  are  zero.  Evaluating  these  determinants  and  setting  them  equal 
to  zero  we  are  led  to 

-  o2pY(y-l)  y2k(x2  +  y2)  +  o[2p(pY-1-l)iB(yriy5  +  xnx?)y2k 

(2.16)  -  kynyk(Y-l)pY(*2  +  y2)]  +  B2  ^-4-  nk  1^-4-  nk(Y-l)pY  2(x2  +  y2) 

*  ft 

-  2yk(pY"1-l)(xr)y5  -  x^y^)  ]  -  iknkB£_|_  pV"1(Y-l)(x2  +  y2)xn  -  0 
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Using  (2.12)  and  the  fact  that  0  is  being  assumed  large  the  above  reduces  to 


*  *  2 


I1  py(y-D  y2k  -  e2(fi-|r-)  n2k  [(y-D  py"2  +  ^(py“1-i)] 

q  p 


or 

,  *  *  2  ,  2k  , 

(2.17)  a2  -  (^-)  e2  (£)  (1-M2) 

Equations  (2.17)  indicate  what  we  would  expect.  In  the  supersonic  (hyperbolic) 
region,  M  >  1,  so  a  is  purely  Imaginary  and  small  perturbations  will  remain  small. 

On  the  other  hand  in  the  subsonic  (elliptic)  region,  M  <  1,  so  a  can  be  positive 
indicating  that  a  small  perturbation  once  introduced  would  grow  as  n  increases. 


10 


NOLTR  70-131 


CHAPTER  III 

THE  NUMERICAL  WORK 

We  are  Interested  here  in  developing  a  numerical  scheme  which  could  be  used 
to  determine  nozzle  contours  in  a  uniform  manner  throughout  the  nozzle.  Our 
procedure  will  be  to  introduce  a  finite  difference  method  to  solve  the  equations 
(2.12)  as  an  initial  value  problem  with  some  Mach  number  or  velocity  array 
prescribed  along  the  centerline. 

There  are  many  ways  in  which  one  could  introduce  finite  differences  to 
approximate  the  differential  equations  we  are  considering.  The  one  used  here  has 
been  chosen  for  its  numerical  simplicity.  We  form  a  rectangular  network  of  points 
(?i»  nm)  with  spacings  A£,  An  and  introduce  the  notation  y,  -  y(£  ,  n  ),  etc.  The 
difference  equations  can  then  be  written 


To  solve  (2.12)  using  (3.1)  we  proceed  as  follows.  We  first  set  xj  o  " 

and  y  ■  0.  The  array  ~  is  either  directly  prescribed  or  obtained  from  a  given 
j  j 

Mach  number  distribution  using  the  nondlmenslonal  form  of  equation  (2.3).  With 


this  information  p.  n  can  be  determined  and  from  this  x  .  and  y,  . .  Using 

JfU  J  »  J- 


these  new  arrays  for  x  and  y, 


is  computed  and  the  procedure  continues 
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Inductively.  Thus,  when  all  the  Information  for  subscript  m  and  less  Is  known 
the  x  and  y  arrays  for  m  +  1  can  be  computed.  Zn  the  program  actually  written 
here,  only  every  other  point  was  calculated  at  each  step.  Thus,  for  example 
(see  Figure  2)  the  point  10  is  calculated  from  1,  2,  and  3,  then  the  points  12 
and  14  would  also  be  calculated  but  not  the  points  11,  13,  and  15.  Similarly, 
going  to  the  next  step,  the  points  10;  11,  and  12  would  be  used  in  calculating 

19.  The  point  21  would  also  be  computed,  but.  not  the  point  20.  Note  that  with 

each  step  in  the  n  direction  1  point  is  lost  from  each  side  of  the  £  array. 

More  will  be  said  about  this  in  Section  IV. 

The  only  deviation  from  this  procedure  occurs  on  the  initial  step  of  the 
process  for  the  axisymmetric  flow  case.  Here  it  is  necessary  to  introduce  the 
approximation 


X 

n 


% 


lx 

3n 


to  circumvent  the  difficulty  arising  from  the  lndetermlnancy  of  along  n  ■  0. 
Using  this  approximation,  we  obtain  for  the  first  step 


(3.2) 


IX 

3n 


1/2 


(since 


1  along  n 


0)  or  in  the  finite  difference  form 


(3.3) 
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The  Iterative  procedure  outlined  above  Is  terminated  when  a  desired  stream¬ 
line,  rij.  Is  reached.  In  the  work  described  here  this  was  chosen  from 
considerations  allowing  comparison  with  the  results  of  others.  As  an  example 
of  how  cine  might  determine  we  consider  the  following  case.  It  was  desired 
to  compute  nozzles  with  nondimensional  test  section  height  equal  to  1.  From  (2.4) 
we  have 


where  pu  Is  constant  in  the  test  section.  Thus,  integrating  from  the  centerline, 
y  ■  0,  to  the  contour  in  the  test  seciton,  y  ■  1,  one  finds 


or  from  (2.11) 


*  * 

p  q 


k+l 


k+1 


- 

k+l 


Thus,  it  follov-T  that 

1 

/  \  <+l 

(3.4) 

This  quantity  can  be  determined  using  e.g.,  [  l]  once  the  test  section  Mach 
number  io  specified. 
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CHAPTER  IV 

RESULTS  AND  CONCLUSIONS 

The  procedure  outlined  in  the  previous  section  was  programmed  and  run  on  an 
IBM  7090.  The  program  consisted  of  a  main  program  which  called  upon  a  subprogram 
to  generate  the  centerline  velocity  array  and  then  carried  out  the  integration  of 
the  equations. 

For  the  two-dimensional  equations  the  cases  run  were  chosen  so  as  to  allow  a 
comparison  with  the  report  of  Baron  [  2  ] •  In  the  examples  considered  there  was 
agreement  in  the  contour  coordinates  to  between  three  and  four  significant  figures. 
These  results  seemed  most  encouraging  especially  when  we  consider  the  crude 
differencing  we  used  and  the  instability  we  know  to  be  present. 

As  we  attempted  to  move  further  upstream  into  the  subsonic  region,  the 
instability  did  become  a  factor  and  wild  oscillations  along  the  calculated  contour 
became  apparent.  This  is  the  behavior  one  would  expect.  If  we  recall  the 
perturbations  we  made  when  we  considered  the  differential  equations,  they  were  of 
the  form  6yQei^^  +  “n  where,  by  (2.17), 


(1  -  M2) 


Thus,  for  fixed  n  we  must  move  sufficiently  upstream  from  the  sonic  line  for  a, 
and  consequently  ean,  to  become  large  enough  to  distort  the  calculations.  We 
point  out,  however,  that  as  we  make  our  centerline  Mach  number  array  pass  through 
M  ■  1  with  steeper  slope,  our  calculations  will  break  down  closer  to  the  sonic 
line. 

In  the  axisymmetric  case  we  compared  our  results  in  the  supersonic  region  with 
those  of  Glowacki  [  6  ),  who  used  the  method  of  characteristics,  and  in  the 
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transonic  region  with  those  of  Hopkins  and  Hill  [7  ] ,  [8  ] ,  who  used  a  series 
method.  For  the  Glowacki  cases  was  chosen  (as  described  in  the  previous 
section)  so  as  to  provide  a  nondimenslonal  test  section  height  of  1.  The 
coordinates  obtained  agreed  to  about  three  significant  firgures.  In  the  tran¬ 
sonic  comparisons  was  chosen  so  as  to  have  the  nondimenslonal  throat  height 
equal  to  1.  These  were  provided  by  Hopkins  and  Hill  who  determined  them  by 
Iteration.  Unfortunately,  the  n^'s  were  so  large  here  as  to  cause  a  divergent, 
oscillatory  condition  almost  everywhere.  In  an  attempt  to  be  able  to  make  some 
type  of  comparison  a  very  coarse  grid  was  employed  (AC,  An  *  0.1  rather  than 
^  0.01  or  0.001  as  in  the  previously  discussed  cases)  and  a  disagreement  in  the 
results  in  the  third  significant  figure  was  observed.  This  again  was  considered 
quite  good  considering  how  coarse  the  grid  was. 

The  two  primary  difficulties  with  the  method  are  the  instability  in  the 
subsonic  region  and  the  loss  of  the  end  points  at  each  step  (oee  discussion  in 
Chapter  III).  The  former  prohibits  us  from  obtaining  the  entire  subsonic  region 
and  always  leaves  us  suspicious  as  to  how  accurate  are  the  values  we  obtain. 

The  latter  difficulty  is  particularly  Important  in  the  design  of  short  nozzles. 

In  this  case  we  may  not  have  many  mesh  points  to  work  with  and  thus,  by  the  time 
we  reach  the  contour,  much  of  the  region  of  Interest  is  lost.  If  we  attempt  to 

f 

decrease  the  mesh  size  to  provide  us  with  additional  points,  we  must  take  more 
steps  in  the  n  direction  and  thus  are  again  in  the  same  situation. 

In  the  introduction  we  mentioned  that  an  alternate  method  is  desirable  to 
avoid  patching  and  for  use  in  short  nozzles.  The  above  results  indicate 
that  such  a  need  still  exists.  Short  nozzles  and  nozzles  with  steep  Mach  number 
gradients  do  not  yield  themselves  nicely  to  the  method  described  herein.  For 
long  nozzles  where  the  Mach  number  gradient  is  not  large  this  method  does  seem  to 

D 
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provide  a  method  of  getting  the  transonic  and  supersonic  regions  In  a  uniform 
way.  Some  pa  :chlng  must  still  be  done,  however,  to  get  the  contour  further 
upstream  In  the  subsonic  region. 
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LIST  OF  SYMBOLS 

velocity  of  sound 

dimension  parameter  (-0  for  two-dimensional;  <*1  for  axiaymmetrlc) 
normalizing  length 
Mach  number  (aq/a) 
pressure 

velocity  vector,  velocity  magnitude 

velocity  components  in  x  and  y  directions,  respectively 

ration  of  specific  heats  (*1.4  for  air) 

streamline  parameter 

potential  line  parameter 

density 

potential  function 
stream  function 
stagnation  condition 
design  condition 
throat  condition 
centerline  condition 
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